ON THE IMPLEMENTATION OF EXPONENTIAL METHODS FOR 
SEMILINEAR PARABOLIC EQUATIONS 
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Abstract. The time integration of semilinear parabolic problems by exponential methods of 
different kinds is considered. A new algorithm for the implementation of these methods is proposed. 
The algorithm evaluates the operators required by the exponential methods by means of a quadrature 
formula that converges like 0(e~'^^/ with K the number of quadrature nodes. The algorithm 
allows also the evaluation of the associated scalar mappings and in this case the quadrature converges 
like 0{e~'^^). The technique is based on the numerical inversion of sectorial Laplace transforms. 
Several numerical illustrations are provided to test the algorithm. 
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1. Introduction. The good numerical results obtained from the application of 
exponential methods to the time integration of stiff semilinear problems, have moti- 
vated much interest on these kind of methods during the last years, see for instance 
[U [51 [31 [HI [ini E] • The problems under consideration can be written in the abstract 
format 

u'{t) = Au{t) + f{t, u{t)), u{0) = uo, < i < T, (1.1) 

where A is a linear operator representing the highest order differential terms and / 
is a lower order nonlinear operator. The solution to the initial value problem p.ip is 
then given by the variation of constants formula and most of the exponential methods 
considered in the literature are constructed from this representation of the solution. 

Let us consider for instance the family of multistep exponential methods devel- 
oped in [2] and briefly reviewed in Section \2A[ Given a stepsize h > 0, n > and 
approximations Un+j ~ u(tn+j), tn+j = {n + j)h, < j < fc — 1, the fc-step method 
approximates the solution u of (jl.H) at tn+k — {n + k)h hy 

k~l 

Un+k = (f>t){k,hA)un + h^(j)j+i{k, hA)A^fn, (1.2) 

3=0 

where /„ = /(t„,M„), A denotes the standard forward difference operator, and for 
A e C and fc > 1, 

Mk,X)^e''\ 0, (fc,A) = ^ e^'^-^)^ ^. ^ ^ ^ da, 1 < j < k. (1.3) 

As we can see in (jl.2p . these methods require the evaluation of 4>j{k, hA), < j < fc, 
for (j)j{k, A) defined in p.3p . This is in fact the main difficulty in the implementation 
of the methods in p.2p and, in general, of exponential methods, since they typically 
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require the evaluation of vector- valued mappings 4){hA), with h the time step in the 
discretization and either 



^(A) = e"'^, A G C, (1.4) 



I e("-")V'^)da, (1.5) 



or 





with m an integer and p{<y) a polynomial. The values of m and p in p.Sp depend on 
the method. For instance, in the case of the methods in (|1.2p . it is clear from (|1.3p 
that m = k, the number of steps of the method, and 

/N f cr \ a{<T-l)...{cr-j + 2) 

For the explicit exponential Runge-Kutta methods constructed in |10j it is m = 1 
and 

?'(^) = (T^' ^ ^ 1' 

as we can see in Section 12.21 (|2.7I) . 

In the present paper we propose a way to evaluate the operators cjj^hA) in (|1.4p 
and (jl.Sp when A in (jl.ip is the infinitesimal generator of an analytic semigroup in a 
Banach space X. Thus, we will assume that A : D{A) C X ^ X is sectorial, i.e., A is 
a densely defined and closed linear operator on X and there exist constants M > 0, 
7 G R, and an angle < S < ^, such that the resolvent fulfils 

\\(zI-A)-'\\<y^, for |arg(z-7)| <^-<5. (1.6) 

|z-7| 

As a side product we also obtain an accurate way to evaluate the mappings (t>{X) 
in (jl.Sp at the scalar level, which is itself a well-known problem in numerical analysis. 
This is exemplified in ^llj with the mapping 

VW^'-^, (1.7) 

which is required for instance by the exponential Runge-Kutta methods of [10, . The 
evaluation of ip for small A by using formula (|1.7p suffers from cancellation error. 
On the other hand, the use of a truncated Taylor expansion only works well for small 
enough A. This implies that for some intermediate values of A it is not very much clear 
how to choose the proper formula and moreover both of them could lose accuracy. For 
A inside a sector | arg(A — 7)] > tt — (5, these difficulties can be overcome by writing 
if in the format p.Sp . with m — I and p(cr) = 1. By doing so, we will be able to 
evaluate <^(A), independently of the size of A, by using essentially the same technique 
developed in principle to evaluate the vector- valued mapping ip(hA). 

In the recent literature, several alternatives have been proposed to evaluate the 
required operators (jj^hA) or alternatively their action on given vector, based on a 
Krylov approach [71 [Sj j on Fade approximants [T] or on a suitable contour integral 
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representation of the mappings by means of the Cauchy integral formula llj. In par- 
ticular, in [TT], the goal is both the evaluation of the mappings 4>{hX) at the scalar 
level, assuming that a diagonalization of A is available, but also at the operator level, 
since it also allows the evaluation of the operators 4>{hA) working with the full matrix 
A. However, despite the good computational results reported in [11], the way of se- 
lecting the parameters involved in the quadrature formulas is not very much clear and 
they depend very strongly on the equation considered and the spatial discretization 
parameters. The algorithm we propose is derived by using Laplace transformation 
formulas and is finally based on a suitable contour integral representation of the map- 
pings (j){hA)^ too. However, the quadrature formulas we use borrow their parameters 
from the method in |14| for the numerical inversion of sectorial Laplace transform, 
where a rigorous analysis of the error is performed together with an optimization 
process to choose the different parameters involved in the approximation. 

In order to apply the quadrature formulas developed in [14j , we derive a represen- 
tation of the operators in (|1.4p and (jl.Sp as the inverses of suitable Laplace transforms 
$(z, hA) at certain values of the original variable a. These Laplace transforms have 
all the form 

^{z,hA)^R{z){zI~hA)-^, (1.8) 

with R{z) a scalar rational mapping of z G C. Due to (jl.6p . the mappings $(2, hA) 
turn out to be sectorial in the variable z, i.e., there exist constants 7 G M and M > 0, 
possibly different from the constants in (II. 6p . such that 



$(z, hA) is analytic for z in the sector | arg(z — 7)] < it ~ 5 and there 

M (1.9) 

\\^(z,hA)\\<- — , for some 1. 

In this way, we reduce the problem of computing (j)[hA) to the inversion of a secto- 
rial Laplace transform ^{z,hA) of the form (|1.8|) . We then use the method for the 
numerical inversion of sectorial Laplace transforms developed in [14j and reviewed in 
Section [31 The inversion method consists on a quadrature formula to discretize the 
inversion formula (|3.3p , so that we finally approximate 

K K 

cj){hA)^ Wie^''Hze,hA) = me'^''Rize){zd~hA)-\ (1.10) 
e=-K e=-K 



with the quadrature weights we and nodes ze given in (j3.6p . The convergence results 
in [M] assure an error estimate in the approximation ()1.10p at least like 0(6""^^/^"^). 
Further, by following this approach, our selection of parameters in the implementation 
of exponential methods depends only on d in ()1.6p . being independent of h and M. In 
the particular case that we want to evaluate at a scalar A in the sector | arg(A — 7) | > 
TT — 5, the approximation becomes simply 

0(A) « V z«,e™^^<f(z,,A)= V wee"''^^^ (1.11) 

ze — ^ 

and the convergence rate improves to an 0{e~'^^). 

The approximation in (jl.lOp allows the computation of all the operators (f>{hA) 
required by an exponential method before the time-stepping begins, so that only the 
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products matrix X vector need to be implemented at every time step. However, for- 
mula (|1.10p requires the computation and storage of all the inverses {zgl — hA)~^, 
£ = —K, . . . , K. Even if A is a sparse matrix and these inversions can be performed 
efhciently, the storage of the resulting full matrices and the subsequent prod- 

ucts matrix X vector can become prohibitive for large problems. In this situation, 
(|1.10p should be combined with a data sparse procedure to approximate the resolvent 
operators, as it is proposed in [HIS]. 

Another way to avoid the computation and storage of all the resolvents in (|1.10p 
could be the application of the formula to evaluate the action (j){hA)v on a given vector 
V, instead of the full operator (j){hA). In this case, it is not necessary the computation 
of the full inverses {zgl — hA)~^, but the resolution of the linear systems 

(zil — hA)x = V. 

However, solving all these linear systems for all the quadrature nodes Z£ in (|1.10p at 
every time step can become quite expensive and, in this sense, the Krylov approach 
developed in [7] appears as a more competitive alternative. In this situation, only for- 
mula p. lip can be useful, in order to evaluate the mappings (j){hX) at the eigenvalues 
of the Hessenberg matrices involved in the Krylov approximation. 

By using (jl.lOp we are in fact computing an approximation to the numerical 
solution of (|l.ip provided by an exponential method. Thus, the global error after 
applying (jl.lOp to the time integration of (|l.ip can be split into the error in the time 
integration by the "pure" exponential method and the deviation from the numerical 
solution introduced by the approximation (jl.lOp of the operators (f){hA). The error in 
the time integration for the exponential integrators considered in the present paper is 
analyzed in |2] and [1^ (see Section [2]), and the quadrature error is analyzed in [14] 
(see Section|31 Theorem l3.1l) . In order to visualize the effect of this approximation, we 
show in Section [6] the performance of our implementation for several problems with 
known exact solution and moderate size after the spatial discretization. In the error 
plots provided we can observe that the error coincides with the expected error for the 
exact exponential integrators up to high accuracy for quite moderate values of K in 
(|1.10p . i.e., the error induced by the quadrature (jl.lOp is negligible compared to the 
error in the time integration. 

Finally we notice that the matrix exponential e*'^ and also certain rational appro- 
ximations to it originating from Runge-Kutta schemes have already been successfully 
approximated by using this approach [12l HH [14] . 

The paper is organized as follows. In Section [2| we consider the class of explicit 
multistep exponential methods proposed in [5] and the explicit exponential Runge- 
Kutta methods in [10]. Section[3]is a review of the method for the numerical inversion 
of sectorial Laplace transforms presented in [M] . In Section [4] we deduce a represen- 
tation for the operators required in the implementation of these integrators in terms 
of suitable Laplace transforms and apply the inversion method to the implementation 
of exponential methods. In Section [5| we consider with some detail the evaluation of 
the associated mappings at the scalar level and present some numerical results. We 
finally test our algorithm with several examples in Section [6l where we implemented 
(jl.lOp by using the full matrices. 

2. Exponential methods. In this section we review some of the exponential 
methods in the recent literature. In this way, we consider the class of explicit multistep 
exponential integrators developed in [2] and the exponential Runge-Kutta methods 
of [TU], which are explicit as well. 
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2.1. Multistep exponential methods. In a class of explicit exponential 
methods is constructed for the time integration of problems of the form Ijl.ip with 
A : D{A) C X ^ X the infinitesimal generator of a Co-semigroup e*^, i > 0, of linear 
and bounded operators in a Banach space X . 

As we already stated in the Introduction, we will restrict ourselves to the case of 
A in p.ip sectorial. Then, for < a < 1, the fractional powers (w — A)" are well 
defined for w > 7 in (|1.6p . and Xa — D{{uj — A)") endowed with the graph norm 
II ■ II Q is a Banach space [3]. In this situation, the nonlinearity / in (|l.ip is assumed 
to be defined on [0, T] x Xa ^ X, for some < a < 1, and to be locally Lipschitz in 
a strip along the exact solution. Thus, there exists L{R,T) > such that 

\\fit,v)-f(t,m<L\\v-ao., V,^^Xa, 0<t<T, (2.1) 

for max (II77 - u{t)\\a, U " w(t)||a) < R- 

For fc > 1, the fc-step method is constructed from the variation of constants 
formula in the interval [t„, i„+fc] as follows. Taking a stepsize h = T/N, N > k, and 
the corresponding time levels t„ = nh, < n < N, the solution u of (|l.ip at tn+k is 
given by 

u{tn+k) - e''^'^u{tn) + h / e('=-")''^/(i„ + ah, u{tn + ah)) da. (2.2) 

Given approximations Un+j ~ uit^+j), < j < /c — 1, the approximation w 
u(tn+k) is obtained after replacing / in (|2.2p by the Lagrange interpolation polynomial 
of degree fc — 1, Pn,,k-i through the points {(t„+j, f(tn+j, Un+j))}j=o ^'^"-^ integrating 
exactly. Writing 

k-i , . 

F„,fc_i(t„ +oh)=Y, \ I AVn, (2.3) 

with /m = fitrm i*m), < m < iV — 1, and A the standard forward difference operator, 
the approximation is computed from 

k-l 

(fc, hA)un + hJ2 <Pj+i{k, hA)A^U, (2.4) 

3=0 

where, for A G C, fc > 1, and < j < fc, the mappings 0j(fc, A) are given in (|1.3p . The 
methods defined in (|2.4p are explicit and, as we already mentioned in the Introduction, 
they require the evaluation of 0j(fc, hA). 

In [21 Theorem 1] it is shown that if the nonlinearity f{t,u{t)) belongs to the 
space C^dO, T],X) and the starting values uq, . . . , Uk-i G Xa satisfy 

\\uitj)-Uj\\a<Coh'', 0<j<fc-l, 

the method defined in (|2.4p exhibits full order fc, i.e., there exists C > such that 

\Ht^)-u„\\a<Ch^\\f<-''^\\oo, 0<n<N. (2.5) 

The constant C depends on k,a,T, L{R,T) in (|2.ip . and 7,M in p.6p . but it is 
independent of ft, and /. 
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2.2. Exponential Runge— Kutta methods. Explicit exponential Runge-Ku- 
tta methods are presented in [lOj for the time integration of semilinear parabolic 
problems. For h = T/N, N > 1, and 1 < i < s, the approximations u„ to u{tn), with 
t„ = nh, are given by 

Um = e^'^^Un + hY^ Qrj {hA)f{tn + Cj/l, Unj), 

(2.6) 

Un+l = e^^Un + hy^bi{hA)f{tn + Cjh, Uni), 

1=1 

with ci = {Uni = Un). The coefficients bi{X) and ay(A) are linear combinations of 

ifkiX) and (fkiciX) with 

Vk{\)= I e^'-'^ ^,""" ' rf^, AeC, A;>1, <> 0. (2.7) 

Setting V'o(-^) = e^, we see that the implementation of (j2.6p requires the evaluation 
of ipk{hA) and (pk(cihA), for 1 < ^ < s and several values oi k > 0. 

As in 2J , the nonlinearity / in is assumed to satisfy a local Lipschitz condi- 
tion like (|2.ip and the solution it : [0, T] — * Xq, and / : [0, T] x Xq, ^ X are assumed 
to be sufficiently smooth so that the composition 

/* : [0, T] X X„ -> X : i /*(i) = /(t, u{t)) 

is a sufficiently often differentiable mapping, too. Under these assumptions, stiff order 
conditions are derived and exponential Runge-Kutta methods of the form (j2.6p are 
constructed up to order four in [10^. 

3. The numerical inversion of sectorial Laplace transforms. In this sec- 
tion we review the numerical inversion method for sectorial Laplace transforms pre- 
sented in [T3], which is a further development of |13j. 

For a locally integrable mapping / : (0, c») —> X, bounded by 

\\f{t)\\<Ce-^e'^\ for some 7 e R, t^ > 0, (3.1) 
we denote its Laplace transform 

poo 

F{z) ^ C[f]{z) = e-''f{t)dt, Rez>7. (3.2) 
Jo 

When F satisfies (|1.9p . the method in jTl] allows to approximate the values of / from 
few evaluations of F. This is achieved by means of a suitable quadrature rule to 
discretize the inversion formula 

fit)^^ f e^'F{z)dz, (3.3) 
2711 Jy 

where F is a contour in the complex plane, running from —ioo to ioo and laying in 
the analyticity region of F. Due to (|1.9p . F can be taken so that it begins and ends 
in the half plane Rez < 0. Following [14], in (|3.3p we choose F as the left branch of 
a hyperbola parameterized by 



M ^ F : X 1-^ T{x) = fi{l — sin(Q! -I- ix)) + 7, 
6 



(3.4) 



Fig. 3.1. Action of T in 13AI on the real axis 



where > is a scale parameter, 7 is the shift in p.9p . and < a < ^ — S. Thus, 
r is the left branch of the hyperbola with center at (/^, 0), foci at (0,0), (2/Lt, 0), and 
with asymptotes forming angles ±(7r/2 + a) with the real axis, so that F remains in 
the sector of analyticity of F, \ a.rg{z — j)\ < tt — 5. In Figure [5?T] we show the action 
of the conformal mapping T on the real axis. 

After parameterizing (|3.3p . the function / is approximated by applying the trun- 
cated trapezoidal rule to the resulting integral along the real axis, i.e., 

K 

fit) = ^(^) ~ E ^(^^)' (3.5) 

-^^ e=-K 

with quadrature weights wi and nodes zi given by 

we ^ --^T'Ut) , ze^Tier), -K < £ < K, (3.6) 

and T > a suitable step length parameter. We notice that the minus sign in the for- 
mula for the weights comes from setting the proper orientation in the parametrization 
of F. In case of symmetry, the sum in p.5p can be halved to 

/(0« Re (5]u;|e*-^i^(z,)) , (3.7) 
\i=o ) 

with Wq = wq and — 2u'^, £ > 1. The good behavior of the quadrature formula 
p.5p is due to the good properties of the trapezoidal rule when the integrand can be 
analytically extended to a horizontal strip around the real axis [T7|, [T8] . 

During the last years, different choices of contours F and parameterizations have 
been studied for the numerical inversion of sectorial Laplace transforms. Apart from 
the approach in [13l [M] , which is the one we follow, the choice of a hyperbola has 
been studied in [H O [TH [TBI EI] ■ The choice of F as a parabola has been considered 
recently in and we refer also to Talbot's method [121 HO] for another kind of 

integration contour F, with horizontal asymptotes as \z\ — > 00. 
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The error analysis for (|3.5p shows that the role played by the round-ofF errors 
is very important. Several ways to avoid error propagation are proposed in [14] . 
depending on the available information about the errors in the evaluations of F{z() and 
the elementary mappings involved in (j3.5|) . Our algorithm for exponential methods 
uses (jl.lOp to approximate the required operators (j){hA) in \1A\ and (|1.5p . i.e., (|3.5p 
with Laplace transforms of the form (jl.Sp . Thus, the required Laplace transforms 
will be in practice evaluated by performing the inversion of few matrices of the form 
zl — hA, with A being normally a discrete version of the operator in (jl.ip and h 
the time step length. Taking into account that this is likely to be accomplished by 
means of some auxiliary routine, we propose to use in this context the least accurate 
version of the method in [T3], where the error propagation is avoided without using 
any information about the errors in the evaluations of (|3.5p . 

The following result provides an error bound for the proposed version of the 
inversion method like 0(e~'^^/'"^), with K the number of quadrature nodes. 

Theorem 3.1. JJ^ Assume that the Laplace transform F[z) satisfies the sectorial 
condition (|1.9p and let a and d he such that 

Q<a~d<a + d<^-5. (3.8) 

For to > 0, A > 1 and K > I, we select the parameters 

with a{K) = arccosh {KK / sma). 

Then, the error E^lt) in the approximation (j3.5p to f{t) with quadrature weights 
and nodes in (j3.6p is bounded by 

/ -2TrdK/a{K) \ 

||£;,(OII <MnQe-'^M-)t-i [e + ^-—^^^.^j , (3.10) 

uniformly for t G [to, Ato], where M and v are the constants in ()1.9p . e is the precision 
in the evaluations of the Laplace transform F and the elementary operations in (j3.5p . 




1 + sin(Q; + d) 
Tr y (1 - sin(a + d))2''-i 

and 

Q = max{4L(Ato sin(Q! — d)), t + L{Xto sin a)}, 

withL{x) = l-ln(l-e-^). 

In case we have some reliable information about the errors in the computation 
of the matrices {zl — hA)^^, an e-depending selection of r and A improves the error 
bound ()3.10p to 0{e^''^), i.e., a pure exponentially decaying bound in K. In this 
situation, given e, K, and a, d fulfilling p.Sp . the parameters t and ii are given by 

a{e*) 27rdK{l-e*) 

where, for 6 G (0, 1), a{d) is the mapping 

a{9) = arccosh (- ) , (3.12) 

V 1 — sm a / 



and 



9*= min fee2-''^(i-^)/<^) + e-2-'^^«M«)V (3.13) 
ee(oa) V / 

Given K > I, the expression to be minimize in (|3.13p represents the main part of the 
error bound obtained in 14J for any fixed 9 S (0,1). By choosing for every K the 
optimal value 9* in (|3.13p . an error bound like 0{e~'^^) is attained. We notice that 
the error bound stated in Theorem l3.1l is obtained for 9* — 1 ~ 1/K in (|3.1ip . instead 
of 9* in (|3.13p . We refer to [14 for details. There, the case < 1 in p.9p is also 
studied. 

4. Evaluation of the vector- valued mappings. In this section we apply some 
Laplace transformation formulas to obtain a suitable representation of the operators 
<j)j{k,hA), ipj{hA), and ipj{cihA) required in ([2^ and ([2^ . 

Let us denote the Laplace transform of a mapping f{a) by F{z) = C[f]{z), and 
the inverse Laplace transform by f{a) = C^^[F]{<j). 

4.1. Evaluation of the mappings required by the multistep methods. 

For in with 1 < j < k, it holds 

0,(fc,A)=^%('=--)^ ( ^ ) da = C-'[C[fo{-,X)]xC[.fMk), 



where, for cr > 0, 



/o(<T,A)=e-^ and f,{a) = ( ) . (4.1) 



J 



For every j > 1 and z G C, wc define 



<i>,(z,A) = /:[/o(-,A)](z) X Cmz) = X C[f,]iz). (4.2) 

Then, for every A G C and j > 1, 

(/..(fc. A) =£-![$,(., A)](fc). (4.3) 

For j = 

„fcA — 1 ' 1 



and thus we define 



<I>o(z,A) = ^. (4.4) 

Z — A 



For A scalar, the mappings ^j{z, A), with 1 < j < 4 are given by 
<i>i(z,A) = — <i>2{z,X) ^ 



z(z-A)' ' ' z2(z-A)' 

* ^ \^ 2-z ^ , 3-3z + z2 

2z'^[z — A) 6z^[z — A) 
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(4.5) 



In order to evaluate 4>j{k^hA), < < 4, we propose to use the formulas in 
with hA instead of A and perform the inversion of the Laplace transform to 
approximate the original mappings at a — k. In this way, the Laplace transforms we 
need to invert are: 

<i>oiz,hA) = {zI-hA)-\ 
<Pi{z,hA) = ^{zI-hA)-\ 

<^2{z,hA) = ^{zI~hA)-\ (4.6) 
<p3{z,hA) = ^-l(zI~hA)-\ 

Mz,hA} = ^-f^ + '\ zI^hA)-\ 

Although the formulas in (|4.6p are derived just formally, we notice that they can be 
justified by combining the Cauchy integral formula with the inversion formula for the 
Laplace transform. More precisely, for suitable contours Fi and r2 in the complex 
plane, both laying in the resolvent set oi A, it holds 



(j)i{k,hA) = -L /" (j)Jk,\){\I ~ hA)-UX 
2m J-p^ 



3^ 

T2 



^ f e''i'S>,{^,hA)d^ = C-'[<S>,{;hA)]{k). 



(4.7) 



Due to (ll.6|l . all the Laplace transforms ^j{k,hA) in (14. 6|) are sectorial, since 
they satisfy (II. 9p with 7* = max{0,7} and 6* — 6, for 7 and S in (|1.6p . We notice 
that, for all j, the resulting bounds in p.9p are independent of h. 

Thus, we can compute the operators 4>j{k,hA), < j < fc, by using the method 
described in Section [3] to compute the inverse Laplace transforms of the mappings 
$j(z, hA) in (|4.2p . We notice that the inverse Laplace transforms need to be approx- 
imated only at the fix value a = k, which is specially favorable for the application of 
the inversion method (see the bound in Theorem 13. ip . Then, we set A = 1, Iq = k, 
and select the parameters r and /i following Theorem 13. II The selection of a and d 
is more heuristic and a good choice is a ~ ^ S) and d slighly smaller than a. 
For example, if 5 = in (|1.9p . good values are around a = 0.7 and d = 0.6. Next, 
we compute the quadrature weights wg and nodes z^ in p.6p and approximate the 
operators in (|2.4p by 



K 



(j)j{k,hA)K wie''^'<^j{ze,hA). (4.: 



e=~K 



The sum in (14. 8p can be halved in case of symmetry like in ()3.7p . 

As we already mentioned in the Introduction, the computation of all the required 
operators (j)j{k, hA), < j < k, can be carried out before the time stepping of 
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the exponential method begins. Thus, if we use the method of Unes and apply the 
exponential method to some spatial discretization of only the matrix- vector 

products in (|2.4p need to be computed at every time step. 

4.2. Evaluation of the mappings required by the Runge Kutta meth- 
ods. For ipj in (|2.7[) . j >1 and t > 0, we have 

where, for ct > and A e C, 

go{a,\) = e'^^ and gj{a) = j—^y- (4-9) 
For every j > 1 and z G C, we define 

^,{z,\) = C[g^{;X)]{z) X C[g,]{z) = "^^^^ (4.10) 

and Vl/o(z, A) = (z — A)^^. Then, for every A G C and j > 0, 

^,(A) = £-1 [«-,(•, A)] (1). (4.11) 

The same argument as in ()4.7p justifies the computation of the operators iy9j(/iA) 
and ipj{cihA), j > 0, 2 < Z < s, by performing the inversion of the Laplace transforms 

^jiz,phA)^^{zI-phAy\ j>0, /3 = 1, Q, (4.12) 

to approximate the original mappings at cr = 1. If (|1.6p holds, the Laplace transforms 
in (|4.12p are also the sectorial in the sense of (|1.9p and we can use the inversion 
method of jl4j . 

As in the case of the methods in (|2.4p , the computation of all the required oper- 
ators ipj{hA) and ipj{cihA), j > 0, can be carried out before the time stepping of the 
exponential method begins. 

Remark 1. In general, we can always evaluate a mapping (j){hA) of the form of 
(jl.Sp by using the numerical inversion of the Laplace transform, just by noticing that 
(j){hA) is the inverse Laplace transform at a ~ n of a mapping ^{z, hA) like in (jl.Sp . 

$(z, HA) = R{z){zl - hA)-\ 

with R{z) = C[p]{z), a scalar rational function of z. 

The above Remark implies that our algorithm can be used to implement other 
kinds of exponential methods, different than those in [51 [TD], as long as they require 
the evaluation of mappings of the form of (|1.4p and (|1.5p . 

5. Evaluation of the scalar mappings. As we already mentioned in the In- 
troduction, we can also apply the inversion of the Laplace transform to evaluate with 
accuracy the scalar mappings 4>{\) in (|1.5p . In this section we consider with some 
detail the evaluation of the mappings 

gj{m,X)::^ e^^'-^^^a^-^ da, j > 1, m e N, (5.1) 
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by means of the quadrature formula p. lip . The result provided by (|l.lip does not 
depend on the size of A, but the formula is only useful in principle for values of A 
inside a sector of the form | arg(A — 7)] > n ~ 5. However, using that 

e'"VK-A)=5iKA), AeC, (5.2) 

and 

g^.^,(™,A) = M!!l^l^, J>1, meN, (5.3) 
it is easy to see by induction that, for m G N and AeC, 

e"^g,(m, -A) = ^ ( £ 1 1 ) A), j > 1. (5.4) 

1=1 ^ ^ 

Thus, we can compute 

g,(m, -A) = e-"^^C-' [G*(-, A)] (m), (5.5) 

with 

G;(^,A) = -^-^y:( ]e\i-iY-\mzy-', j>i. (5.6) 



which provides a stable formula to approximate gj{m, — A) for A inside a proper sector 
I arg(A — 7)1 > TT — S and moderate size. 

In Table [01 we show the error obtained in the evaluation of (p{X) — gi{l,X) 
in (|1.7p for different values of A in the interval [—1,1]. For A < 0, we applied the 
inversion formula (j3.5p with t — I and 

F(z) = Gi(z,A) = — (5.7) 
z(z — A) 



which, for these values of A, fulfils ()1.9I) with 6 = 0,^ — 0, and i> — 2. We assumed 
that the evaluations of Gi can be carried out in MATLAB up to machine accuracy 
and thus we set e = 2.2204 x 10~^^. Then, we computed the quadrature weights and 
nodes in following ((3TT|) - ((3T3)) with A = 1. Setting a = 0.7 and d = 0.6, we 
obtained = 0.693, for K = 15, and 9 = 0.793, for K = 25. In Table O we can 
see that ii' = 25 is enough to attain almost the machine accuracy of MATLAB in the 
evaluations of ip{X). For positive values of A, we used (|5.2p with m = 1. 

6. Numerical illustrations. In this section we test our algorithm by consider- 
ing the same examples as in 2 and pUj . 

6.1. Example for the multistep exponential methods. Our first example 
is the problem considered in [2] 

ut{x,t) ^ Uxx{x,t) + (^j u{s,t)ds^ Ux{x,t) + g{x,t), (6.1) 

for X e [0,1] and t G [0,1], subject to homogeneous Dirichlet boundary conditions 
and with g{x,t) such that the exact solution to (j6.ip is u{x,t) = a;(l — x)e* 
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Table 5.1 

Computation of ^p{X) in 111.71 1 for X € [—1, 1] by using formulas i3.5t and 115.21 1 . We show the 
absolute error obtained in MATLAB with K = 15 and K = 25. 



A < 




K = 1 


5 


K = 2b 


-A 


K = 15 


K = 25 


-1 


1 


5050e- 


12 


1 


3323e- 


15 


1 


1 


5050e- 


12 


3.3307e-15 


-le-1 


1 


5227e- 


12 


3 


2196e- 


15 


le-1 


1 


5227e- 


12 


3.5527e-15 


-le-2 


1 


4243e- 


12 


4 


4409e- 


15 


le-2 


1 


4243e- 


12 


4.6629e-15 


-le-3 


1 


3750e- 


12 


1 


3323e- 


15 


le-3 


1 


3750e- 


12 


1.3323e-15 


-le-4 


1 


3738e- 


12 


1 


7764e- 


15 


le-4 


1 


3738e- 


12 


1.7764e-15 


-le-5 


1 


3747e- 


12 


3 


6637e- 


15 


le-5 


1 


3747e- 


12 


3.7748e-15 


-le-6 


1 


3748e- 


12 


3 


6637e- 


15 


le-6 


1 


3748e- 


12 


3.7748e-15 


-le-7 


1 


3695e- 


12 


1 


9984e- 


15 


le-7 


1 


3695e- 


12 


1.9984e-15 


-le-8 


1 


3717e- 


12 


1 


1102e- 


16 


le-8 


1 


3717e- 


12 


2.2204e-16 


-le-9 


1 


3715e- 


12 


1 


1102e- 


16 


le-9 


1 


3715e- 


12 





-le-10 


1 


3711e- 


12 









le-10 


1 


3711e- 


12 





-le-11 


1 


3711e- 


12 









le-11 


1 


3711e- 


12 





-le-12 


1 


3715e- 


12 


1 


1102e- 


16 


le-12 


1 


3715e- 


12 





-le-13 


1 


3712e- 


12 









le-13 


1 


3712e- 


12 


2.2204e-16 



The spatial discretization of (j6.ip is carried out by using standard finite differences 
with J = 512 spatial nodes, centered for the approximation of Ux- The nonlocal term 
is approximated by means of the composite Simpson's formula. 

To integrate in time the semidiscrete problem we use (|2.4p with k = 1,2,3 and 4, 
so that A is the ( J — 1) x ( J — 1) matrix 

A = J^tridiag ([1,-2,1]). 

We approximate the matrices 0j(fc, hA), < j < k, required in (|2.4p by applying the 
quadrature rule ()4.8|) . To avoid an extra source of error, the initial values wi, . . . , Uk-i 
are computed from the exact solution. In a less academic example, these values can be 
computed by means of a one-step method of sufficiently high order or by performing 
the fix point iteration proposed in [2] . 

In Figure 16.11 we show the error versus the stepsize at t — I, measured in a 
discrete version of the norm || • ||i/2i for K = 25 and X = 35 in (|4.8p . We see that for 
if = 35 the full precision is achieved for all the methods implemented; cf. Section 
6]. In Figure [Ol we also show lines of slope 1, 2, 3 and 4, to visualize the order of 
convergence. In fact, we can observe that the order of convergence for this example 
is slightly higher than the one predicted in approximately 2.15, 3.15 and 4.15 for 
A; = 2, 3 and 4, respectively, instead of sharp order k. A further study of this behavior 
is beyond the scope of this paper. 

6.2. Examples for the exponential Runge Kutta methods. For the fol- 
lowing two examples we consider the problems and some of the integrators proposed 
in [10 . Following the notation in [TU], in the Butcher tableaus below we use the 
abbreviations 

ipi = (Pi{hA), and ipi^j = (ptj{hA) = (p,{cjhA), 2 < j < s. (6.2) 
Thus, our second example is 

ut{x,t) ^ Uxx{x,t) + — — - — —+g{x,t), (6.3) 

1 + u[x, ty 
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Fig. 6.1. Error of exponential multistep methods 112. 4| I applied to 116.11 1 . for k = 1, 2, 3, and 4. 
Left: With K = 25 quadrature nodes on the hyperbolas, Right: With K = 35. 



for X G [0,1] and t G [0,1], subject to homogeneous Dirichlet boundary conditions 
and with g{x,t) such that the exact solution to (|6.3p is again u{x,t) = x{l — a;)e*. 

We discretize this problem in space by standard finite differences with J = 200 
grid points. For the time integration of the semidiscrete problem, we implemented 
(12.61) with s = 1, the second-order method 








1 




2 






VI 



(6.4) 



the third-order method 










1 

3 


5<^1,2 




2 
3 


§<<5l,3 - |<^2,3 


|<^2,3 




fl - |V2 


1^2 



(6.5) 



and the fourth-order one 



with 












1 








2 


5'y3l,2 






1 
2 


- V2,3 


</'2,3 




1 


(y9i,4 - 2(y92,4 


</?2,4 ¥^2,4 




1 
2 


^</'i,5 — 2a5^2 — as, 4 


15,2 ^5,2 15,4 








-<y92 + 4(^3 4932 — i 






1 


1 1 






05,2 = 2'/'2,5 


- '/'3,4 + ^"^2,4 - -^^3,5 





(6.6) 



and 



05,4 



1 



r</52,5 - 05,2- 
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For the implementation of (|6.4p we need to invert four different Laplace trans- 
forms of the form of (|4.12p . to approximate ^po (f^) T^oi^A), ipi (f^), and ipi{hA). 
The implementation of both (j6.5p and (|6.6p requires the inversion of eight Laplace 
transforms. 

In Figure 16.21 we show the error at t = 1 versus the stepsize, measured in the 
maximum norm. The expected order of convergence for this example is k for the k- 
order method. In order to check our algorithm, we added lines with the corresponding 
slopes in Figure 16.21 We can see that also for this kind of methods we attain full 
precision for if = 35 in (|4.8p . 

Finally we consider the second example in [TU] 

ut{x,t) = Uxx{x,t) + u{s,t) ds + g{x,t), (6.7) 
Jo 

for X e [0, 1] and t £ [0, 1], subject also to homogeneous Dirichlet boundary conditions 
and with g{x,t) such that the exact solution to (|6.7p is u{x,t) — x{l — x)e*. 

We discretize this problem in space as in the previous example (|6.3|) , and use the 
composite Simpson's rule for the approximation of the nonlocal term. For the time 
integration, we use (|2.6p with 5 = 1, (|6.4p . and (|6.5p . In Figure [^751 we show the error 
at t = 1, measured in a discrete version of the norm || • II1/2; cf. flUl Section 6]. The 
expected order of convergence for this example is 1 for the first order method, 1.75 for 
()6.4p . and 2.75 for ()6.5p . In order to test our algorithm for the exponential methods, 
we added lines with the corresponding slopes in Figure 16.31 We refer to ^lOj for a 
detailed explanation of the order reduction phenomenon in this example. The less 
stringent accuracy requirements in this case allow us to achieve full precision with 
only if = 25 quadrature nodes in (14. 8p . 

7. Conclusions. In this paper we derived a way to approximate the exponential- 
like operators required for the implementation of different kinds of exponential me- 
thods present in the recent literature for the time integration of scmilinear problems. 
The approach is based on the numerical inversion of the Laplace transform and its 
application is restricted to parabolic problems. When applicable, the proposed algo- 
rithm is shown to be very efficient, both at the scalar level and the operator level, and 
it is quite simple to implement. Apart from the error bound stated in Theorem 13.11 
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Fig. 6.3. Error of Runge-Kutta methods 112.61 1 with s = 1, I I6.4I I . and 116.51 1 . applied to 16. 7D . 
Left: With K = 20 quadrature nodes in Ij4.8|l . Right: With K = 25. 

which is a partial result of those proved in [T^, we tested the algorithm with sev- 
eral academic examples from [2 and [10 and with the evaluation of the prototypical 
mapping ip in p.7p . 

The comparison of our algorithm with other techniques present in the literature 
(see the Introduction) is beyond the scope of the present paper. Some testing for com- 
parison of CPU times and storage requirements with problems of different complexity 
could be the subject of future work. 
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